##### for Figure 4 a,c

rm(list=ls())
setwd("/Users/gin/Documents/onedrive/data/public_database/scRNA/GSE167002/")

library(Seurat)
##### data format ##### -----

### expression matrix
raw.data <- read.csv("GSM5090680_SpleenRC_NI_assay-data.csv", header=T, row.names = 1)

### meta data
metadata <- na.omit(read.csv("GSM5090680_SpleenRC_NI_metadata.csv", header=T))

sce <- CreateSeuratObject(counts = raw.data)

rownames(metadata) = rownames(sce@meta.data)
metadata = metadata[,-1]

sce <- AddMetaData(object = sce, metadata = metadata)

sce <- NormalizeData(sce, normalization.method =  "LogNormalize", scale.factor = 10000)

sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 3000) 

sce <- ScaleData(sce) 

sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce)) 

sce <- FindNeighbors(sce, dims = 1:15)

sce <- FindClusters(sce, resolution = 0.1)
sce <- AddMetaData(object = sce, metadata = metadata)

sce <- RunTSNE(sce)
sce <- RunUMAP(sce,dims = 1:15)
DimPlot(object = sce, reduction = "umap",label = TRUE)
DimPlot(object = sce, reduction = "tsne",group.by= "seurat_clusters",label = TRUE)
DimPlot(object = sce, reduction = "umap",group.by= "seurat_clusters",label = TRUE)
DimPlot(object = sce, reduction = "umap",label = TRUE,pt.size = 1)


FeaturePlot(sce,"Grem1",order=T,pt.size = 1)




##### stackvlnplot ##### -----
library(Seurat)
library(ggplot2)
### Grem1+ and Grem1- ###
g_pos = subset(sce,ident = 0)
g_neg = subset(sce,ident = c(1:6))

meta_gpos = as.data.frame(colnames(g_pos))
colnames(meta_gpos) = "cell"
meta_gpos$group="g_pos"

meta_gneg = as.data.frame(colnames(g_neg))
colnames(meta_gneg) = "cell"
meta_gneg$group="g_neg"

meta = rbind(meta_gpos,meta_gneg)
meta_sce = as.data.frame(colnames(sce))
colnames(meta_sce) = "cell"

meta = merge(meta_sce,meta,by="cell",sort = F)
rownames(meta) = meta[,1]
#meta = meta[,-1]

sce <- AddMetaData(object = sce, metadata = meta[,2], col.name = "group")


### 运行
#genelist = c('Grem1',"Ccl21a",
#             "Tcf21","Wt1","Tlx1","Ly6a",
#             "Dpt","Gsn","Thy1","Il33",
#             "Madcam1","Cxcl13","Igfbp3")
genelist = c('Grem1',
             "Tcf21","Tlx1",
             "Bst1","Ccl19","Ccl21a","Crym","Il33",
             "Madcam1","Cxcl13","Igfbp3")
VlnPlot(sce, features = genelist, pt.size = 0 ,stack=T,group.by = "group") 
FeaturePlot(sce,c("Vim","Dcn"),order=T)
#####

rm(list=ls())
setwd("/Users/gin/Documents/onedrive/data/public_database/scRNA/GSE167002/")

library(Seurat)
##### data format ##### -----

### expression matrix
raw.data <- read.csv("GSM5090680_SpleenRC_NI_assay-data.csv", header=T, row.names = 1)

### meta data
metadata <- na.omit(read.csv("GSM5090680_SpleenRC_NI_metadata.csv", header=T))

sce <- CreateSeuratObject(counts = raw.data)

rownames(metadata) = rownames(sce@meta.data)
metadata = metadata[,-1]

sce <- AddMetaData(object = sce, metadata = metadata)

sce <- NormalizeData(sce, normalization.method =  "LogNormalize", scale.factor = 10000)

sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 3000) 

sce <- ScaleData(sce) 

sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce)) 

sce <- FindNeighbors(sce, dims = 1:15)

sce <- FindClusters(sce, resolution = 0.1)
sce <- AddMetaData(object = sce, metadata = metadata)

sce <- RunTSNE(sce)
sce <- RunUMAP(sce,dims = 1:15)
DimPlot(object = sce, reduction = "umap",label = TRUE)
DimPlot(object = sce, reduction = "tsne",group.by= "seurat_clusters",label = TRUE)
DimPlot(object = sce, reduction = "umap",group.by= "seurat_clusters",label = TRUE)
DimPlot(object = sce, reduction = "umap",label = TRUE,pt.size = 1)


FeaturePlot(sce,"Grem1",order=T,pt.size = 1)




##### stackvlnplot ##### -----
library(Seurat)
library(ggplot2)
### Grem1+ and Grem1- ###
g_pos = subset(sce,ident = 0)
g_neg = subset(sce,ident = c(1:6))

meta_gpos = as.data.frame(colnames(g_pos))
colnames(meta_gpos) = "cell"
meta_gpos$group="g_pos"

meta_gneg = as.data.frame(colnames(g_neg))
colnames(meta_gneg) = "cell"
meta_gneg$group="g_neg"

meta = rbind(meta_gpos,meta_gneg)
meta_sce = as.data.frame(colnames(sce))
colnames(meta_sce) = "cell"

meta = merge(meta_sce,meta,by="cell",sort = F)
rownames(meta) = meta[,1]
#meta = meta[,-1]

sce <- AddMetaData(object = sce, metadata = meta[,2], col.name = "group")


### 运行
#genelist = c('Grem1',"Ccl21a",
#             "Tcf21","Wt1","Tlx1","Ly6a",
#             "Dpt","Gsn","Thy1","Il33",
#             "Madcam1","Cxcl13","Igfbp3")
genelist = c('Grem1',
             "Tcf21","Tlx1",
             "Bst1","Ccl19","Ccl21a","Crym","Il33",
             "Madcam1","Cxcl13","Igfbp3")
VlnPlot(sce, features = genelist, pt.size = 0 ,stack=T,group.by = "group") 

#####

